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Abstract 

We investigate the transport and separation of overdamped particles under the action of a 
uniform external force in a two-dimensional periodic energy landscape. Exact results are obtained 
for the deterministic transport in a square lattice of parabolic, repulsive centers that correspond to 
a piecewise-continuous linear-force model. The trajectories are periodic and commensurate with 
the obstacle lattice and exhibit phase-locking behavior in that the particle moves at the same 
average migration angle for a range of orientation of the external force. The migration angle as a 
function of the orientation of the external force has a Devil's staircase structure. The first transition 
in the migration angle was analyzed in terms of a Poincare map, showing that it corresponds to a 
tangent bifurcation. Numerical results show that the limiting behavior for impenetrable obstacles 
is equivalent to the high Peclet number limit in the case of transport of particles in a periodic 
pattern of solid obstacles. Finally, we show how separation occurs in these systems depending on 
the properties of the particles. 
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I. INTRODUCTION 



The transport of suspended particles in a two-dimensional (2D) periodic periodic lattice 
of obstacles has recently attracted significant attention as a means of separating particle 
species in microfiuidic devices. A driving force induces particles to move through the sys- 
tem. Depending on the properties of the particles and the direction of the driving force 9 
with respect to the lattice orientation, the particles move through the lattice in commensu- 
rate phase-locked trajectories. That is, for a range of driving-force angles single species of 
particles will have an average trajectory in one of the lattice directions [p, g], whose angle 
is given by a = arctan(g/p). Further, that average trajectory angle will remain fixed over 
a range of driving-force angles. When the driving-force angle is outside that range, the 
particle trajectories will jump to a new fixed angle. The relation between the direction of 
motion a and the angle of the driving force 9 depends on the properties of the particles. 
Consequently, the system has the potential to separate different species of particles, with 
the advantage that different particles migrate at different angles {vector chromatography), 
allowing for continuous fractionation [III2]- 

These microfiuidic systems can be broadly classified into two groups, depending on the 
nature of the energy landscape that is experienced by the suspended particles. In one group, 
the particles move through a periodic array of solid obstacles, such as in the separation 
devices based on deterministic hydrodynamics pioneered by Austin's group [SI IH El E]- In 
terms of the corresponding energy landscape, the array of obstacles can be considered as 
a periodic array of hard-core repulsive potentials. In the second class of systems, particles 
move through smooth potential landscapes, with soft interactions between the particles and 
an external field. Examples include the optical fractionation methods pioneered by Grier's 
group, where colloidal particles are transported through an array of holographic optical 
tweezers [3 El El 1101 [II] • 

Phase-locking behavior is common to transport through periodic structures in many 
systems (121 [131 El [13 [E] as well as to non-linear dynamical systems in general [HI [H] . 
In the transport of particles in periodic systems, Lacasta et al. |19] studied a 2D periodic 
arrangement of wells (or traps) in a square lattice by means of numerical simulations of 
the corresponding Langevin equation. Numerical results showed the presence of periodic 
trajectories and lateral migration, in that the particles moved, on average, at an angle a 
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different from the orientation angle of the driving force 9. They also showed the presence of 
phase locking in the a vs. 9 curve, with clear plateaus at large Peclet numbers. In fact, the 
authors identify the observed migration as a deterministic phenomenon [20j. 

In a separate study, Lacasta et al. [21] investigated periodic landscapes that present 
either repulsive (obstacles) or attractive (traps) centers located on a square lattice. In 
both repulsive and attractive cases, they observe similar phase locking behavior, with the 
corresponding plateaus in the a vs. 9 curve becoming evident at large Peclet numbers (Large 
Peclet numbers corresponds to low temperatures in Refs. [IH] and pT]). 

On the other hand, separable potentials, which can be written as a sum of periodic fields 
in each of the two principal directions of a square lattice [2nj, do not exhibit the same 
complex behavior, as clearly shown by Pelton et al. [H]. Experimental work using optical 
tweezers has also demonstrated the presence of periodic trajectories and locked-in states 
[TJ [11] . However, these experiments displayed a different type of locking in the presence of 
Brownian motion, in which the trajectories become commensurate only in a statistical sense 
and only some of the locked-in states are centered on commensurate directions (A hopping 
model for the observed dynamics is proposed by Gopinathan and Grier |9]). 

Conversely, and in agreement with the discussed numerical results at high Peclet num- 
bers, recent simulations and experiments on the motion of non-Brownian spherical particles 
through periodic arrays of obstacles clearly show that deterministic systems exhibit analo- 
gous phase-locking dynamics and periodic trajectories into commensurate lattice directions 

[?]■ 

Phase-locking dynamics has also been observed in more complex transport systems, in- 
cluding the motion of an overdamped particle through a periodic potential under time- 
periodic driving forces [T5| 122] and the dynamics of a particle driven through a deformable 
colloidal lattice |T6] . 

In contrast to a significant number of experiments and simulations, there are few studies 
that present analytical results for the transport of particles in irreducible two-dimensional 
systems. Pelton et al. [H] studied the overdamped motion of a single particle in a periodic 
landscape. In the ID case (linear fringes) they showed the existence of locked trajectories 
for driving forces that are not strong enough to drive the particles over the energy barrier. 
They also showed that the behavior in separable (reducible) 2D potentials is analogous to the 
locking observed in ID systems. Gleeson et al. p3] studied irreducible 2D periodic potentials 
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and derived an iterative method for calculating the average velocity of the particles in inverse 
powers of the external force and particle's diffusivity. 

In this work, we are interested in the deterministic limit (small diffusivity) for relatively 
small forces, i. e., driving forces for which the effect of the potential landscape is not neg- 
ligible. Therefore, we consider both the deterministic evolution as well as the high Peclet 
number limit for finite forces. In particular, we consider a quadratic, continuous model for 
the 2D periodic potential that corresponds to a piecewise-continuous, linear-force model. 
This simple model offers the interesting combination of an irreducible but solvable potential 
that captures the non-trivial transport phenomena exhibited by 2D periodic systems, includ- 
ing directional locking and the universal behavior of dynamical systems near a bifurcation 
point. 

II. TRANSPORT OF COLLOIDAL PARTICLES THROUGH PERIODIC LAND- 
SCAPES 

A. Equation of Motion: High Friction Limit 

The equation of motion for a Brownian particle traversing a periodic force field is the 
Langevin equation [2^, which in the high-friction limit takes the form [25] : 

7^ = F(x) + Fo(x)+^(t), (1) 

where F(x) is the periodic force field, Fq is an external driving force and ^(t) is the Langevin 
force describing the fluctuating force exerted by the fluid on the colloidal particle. The 
friction constant 7 is given by Gir^a, where a is the radius of the colloid and ^ is the 
viscosity of the fluid. The Langevin force is represented by a Gaussian distribution with 
zero mean, (^(t)) = 0, and 6 correlation, (^j(t)^j(s)) = 2'ykT6(t — s)6ij. Here, we shall 
focus on a spatially uniform external force, Fo(x) = Fq and, for convenience, we choose a 
coordinate system with the x axis oriented along the direction of the force, Fo(x) = FqCx 
(see figure [T]) . 
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FIG. 1: (Top) Potential landscape with quadratic, repulsive centers on a square lattice. (Bottom) 
Schematic view of the system, with circles representing the repulsive centers, ei and e2 are the 
principal vector of the square lattice. A unit cell is represented by the dashed line. The Fq is 
oriented along the x axis, and 6 is the angle between the force and ei (in this case is 6 = 25°). 
The solid line shows a locked trajectory that moves in the [2,1] lattice direction, indicated with a 
solid arrow at the bottom of the system, a is the angle between this asymptotic direction and ei: 
a = arctan(l/2) = 26.56°. 
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B. Periodic Potential: Piecewise-linear model 

We consider the case in which the periodic force can be derived from a potential field, 
F(x) = — VK(x). Here, we shall model the two-dimensional, periodic landscape as a 
piecewise-smooth potential that is composed of repulsive centers of size R (obstacles) lo- 
cated in a square lattice with lattice spacing L > R. Specifically, we consider the periodic 
landscape shown in figure [T] which, in the unit cell, is given by. 



where the center of the coordinate system coincides with the center of the cell, r is the 
radial position, = x"^ + y"^, and Fmax gives the magnitude of repulsive centers. Since 
we have chosen to align the x axis with the external force the obstacle lattice will be, in 
general, rotated with respect to the coordinate system. We shall refer to the rotation angle 
between the x axis of the coordinate system and the principal axis of the square lattice ei 
as the forcing angle 9 (see figure [T]). We shall also non-dimensionalize our variables using 
Uc = Fmax/l as the characteristic velocity, Fmax as the characteristic force, and R as the 
characteristic length. The new variables become x' = x.'-f/Fmax', x' = x/i?; r' = r/R. The 
boundary of the repulsive centers is at r' = 1, and we define the relative separation between 
the repulsive centers a.s i = L/R. For simplicity, we do not use the primes to refer to the 
non-dimensional variables in what follows. 

III. DETERMINISTIC TRANSPORT: EXACT SOLUTIONS 

The trajectories in the deterministic limit are obtained from Eqs. [T] by neglecting the 
effect of thermal fluctuations. Thus, outside the quadratic regions, the particles follow a 
straight line that is parallel to the x axis. On the other hand, inside the parabolic regions, 
the equations of motion become. 




r > R, 



r < R 



(2) 



X = 



y = 




r < 1, 



(3) 



with / = Fo/Fmax the ratio of the driving to the repulsive force. 
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We can include the external force in a modified potential field as Vm{x, y) = V{x, y)+x f. 
This modified potential is also an inverted parabola, but with its center shifted to xq = 
(— /, 0). It is clear then that the family of curves perpendicular to the equipotential lines, i. 
e., the particle trajectories inside the circle, are also straight lines, with center at xq. The 



When / > 1 the x component of the velocity is always positive and all the trajectories are 
unbounded. Then, we can define the migration angle a of the particles as the asymptotic 
angle that the trajectory of the particles forms with the principal axis of the lattice ei, i. e., 
a = 9 means that the particles move colinearly with the external force; a = means that 
the particle moves parallel to the principal direction of the square lattice (see figs. [l]and|2]). 

Let us note that we are in the deterministic and high-friction limits, and therefore two 
independent trajectories cannot cross each other in real space. In addition, we consider 
two-dimensional trajectories in the plane, and thus for a given value of the forcing angle 9 
all the trajectories must have the same migration angle a. Therefore, we can obtain a single 
valued a vs. 9 curve by determining the angle of a single trajectory as a function of the 
forcing angle. 

A. Analytical solutions in the unit cell 

In general, a global trajectory can be segmented into a number of successive collisions 
between the colloidal particle and the repulsive centers. In addition, we can classify each 
trajectory passing through a unit cell in terms of its incoming impact parameter 6o which 
we define in the present coordinate system as the y coordinate of the particle when it enters 
the unit cell (see fig. [2^). For 6o > 1 the trajectories are straight lines parallel to the x axis 
and do not interact with the obstacle at the center of the cell. For Kh^ < 1, on the other 
hand, we have three different cases. For impact parameters larger than a maximum value. 



same result can be obtained from direct integration of Eqs. [3j 




it can be shown from Eqs. |3] that the radial component of the velocity is negative outside 
the circle and positive inside it: 




for 



for 



r > 1 



r<l & 0<6o<&; 



(5) 
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FIG. 2: a) Schematic view of a unit cell for a rotation angle 6. bo is the incoming impact parameter 
for the trajectory shown with the solid line. We also show the incoming impact parameter for the 
next collision and the corresponding shift by /?. £ is the ratio between center size R and center-to- 
center distance L. The center for the collision trajectories is (— /, 0). The critical impact parameters 
are be, which corresponds to a collision trajectory ending at (0,±1), and bm, which corresponds 
to a collision trajectory tangent to the parabolic center, b) Schematic view of the three types of 
trajectories: Trajectory A for incoming parameters be < bo < bm', B for < 6o < be, and C for 
< ^0 < 1- ^+ is the maximum value of the incoming parameter for a given forcing direction 
6. The region between and corresponds to the trajectories crossing the top (bottom) of the 
unit cell. 
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In our piecewise approximation, this means that the particle will move around the circle 
with r = 1, until it separates at the top of the circle and then follows a straight line parallel 
to the X axis. A trajectory of this type is shown in fig. |2}d: the incoming particle enters the 
circle at the point Ci with impact parameter bo < —bm, and therefore leaves the circle at 
C2 = (0, —1). The maximum value of the impact parameter, bm, corresponds to the point on 
the circumference for which its tangent passes through the center of the effective potential, 
(— /, 0), as shown in fig. |2^. 

For impact parameters smaller than bm there is a critical value, 

such that, for be < b^ < bm the particle enters the repulsive centers, but the separation from 
the circle still occurs at the top (or bottom for negative values of the impact parameter, 
see below). Such a characteristic trajectory is shown in fig. ^p: the particle enters the 
circle at the point Ai and follows a straight line with center at (— /, 0). Then, the particle 
reaches the point A2 on the circumference of the circle, which belongs to the region above 
bm- Thereafter, the trajectory becomes identical to those described for 60 > bm, with the 
particle separating from the circle at the top (point A^). We will refer to collisions with 
incoming parameter larger than b^ as irreversible in that, independent of the exact value of 
the impact parameter, all the incoming trajectories collapse into a single outgoing trajectory, 
with outgoing impact parameter bf = 1. 

Finally, for impact parameters smaller than the critical value, < 60 < be, the particle 
enters the circle (e. g. point Bi in figure |2}d), moves in a straight line, and leaves the circle on 
the positive side of the x axis (point i?2). The outgoing impact parameter (or y coordinate) 
in this case is given by, 

iP - 1) ^0 



(/2 + i)-2/v/r^ 

B. Poincare Map and Saddle-Point Bifurcation 

In figure |2)d we showed the different types of trajectories that a particle follows depending 
on the impact parameter. For each one of these trajectories, we also show the impact 
parameter that the particle will have in its next collision in a neighboring unit cell (see 
points A\, B[, and C[). If the outgoing impact parameter bf = 6/(&o) is defined as the y 
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coordinate of the point where the particle leaves the cell, then the next incoming parameter 
will he b'^ = bf — P, with /5 = ism(9), as shown in figure [2^. Therefore, for positive values 
of the incoming impact parameter, the impact parameter for the next collision is given by: 

bc<bo<l (8) 
1 < 60 < ^+ 



b' 



(/2+l)-2/^/T^ 

1-/5 
bo- (3 



where = (£/2)(cos(6') + sin(6')) is the maximum possible value of the incoming impact 
parameter (see fig. |2|d). The symmetric conditions apply to negative values of the impact 
parameter, with the only difference being that those trajectories leaving at the bottom of 
the cell come into the next cell from the top, that is: 



b' 



(3 + 2^ 



icos{9) 



<bn < 



(9) 



We can then investigate the global trajectories by studying the above transformation of 
the impact parameter, which is in fact a Poincare map of the impact parameter into itself 
6q = M(6o) [26]. Our objective is to compute the migration angle a as a function of the 
forcing angle 6. We showed before that for a given 6 all the trajectories have the same 
migration angle. Therefore we can investigate single trajectories using the Poincare map 
to determine the asymptotic direction of motion. This is particularly simple in the case of 
periodic trajectories. 

Let us define a periodic trajectory to have period q if the particle goes through q different 
collisions before repeating its motion (Note that all periodic trajectories are commensurate 
with the obstacle lattice) ^7\. This corresponds to q successive collisions with different 
incoming impact parameters before the sequence repeats itself. Periodic trajectories with 
period q correspond to fixed points b^ of the g-times composed map of the impact parameter, 
for which 6q = M^'^^b^). In figure [s^ we present the map of the impact parameter into itself 
for increasing values of 6. We also show the intersection points with the diagonal, which 
correspond to fixed points of period q = 1 (the incoming collision parameter is always the 
same). The presence of fixed points with q = I indicates that the trajectory remains locked 
at a = for sufficiently small forcing angles. In fact, for small forcing angles there is only 
one stable fixed point, located in the region of irreversible collisions (impact parameters 
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FIG. 3: a) Poincare map of the incoming collision parameter into itself for different forcing angles 
e. {I = 2.5; / = 2.0; hm = \/3/2; = 0.6; (ib = 0.5; Oi, = 11.54°.) The map has always the same 
form but is shifted down by a constant amount (3 = ^sin(0). The fixed points corresponding to 
increasing forcing angles are: h\: a fixed point in the region of irreversible collisions; b^'- a fixed 
point that corresponds to successive collisions with bo = be, b^: a fixed point corresponding to 
reversible collisions; b\: a fixed point at the bifurcation angle = Oh- The last map has no fixed 
points corresponding to g = 1. b) A trajectory is shown for a forcing direction (3 = 0.55 > 
The trajectory has a periodicity q = 'd. Only one of the collisions leads to the particle crossing the 
top-bottom boundary of the unit cell and therefore the direction of the trajectory is [8,1]. 
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FIG. 4: Poincare map of the incoming collision parameter into itself and compositions of the map 
of different orders q =1; 2; 3; 5; 7 and 9. The forcing angle is the same as in figm'e [Sja, with 
(3 = 0.55 > (3b. {i = 2.5; / = 2.0; bm = V3/2; be = 0.6; = 0.5; Ob = 11.54°.) It is clear that only 
q = 9 has fixed points, as shown in figure |3jo by following an individual trajectory. 

between be and R). The fixed point corresponds to 6q = 1 — /? {b\ in the figure). The other 
fixed point is unstable, since the local slope of the map is greater than one [TTj. As the 
forcing angle increases, the map is shifted down by an increasing amount, (3 = isin{6), and 
the fixed point moves to the left in the map, corresponding to smaller impact parameters. 
Eventually, the fixed point reaches the critical impact parameter 6q = be, as indicated in the 
figure (62)- For larger forcing angles, the stable fixed point is given by (63 in the figure), 

KiP) = -^ + ^^2p-2 + pp^ + iP - 1) v^W^. (10) 

The corresponding collision not only penetrates the parabolic regions but is also no longer 
irreversible. However, these trajectories are still locked into the a = overall motion. 
Finally, as the forcing angle increases, the map goes through a tangent bifurcation for 
Pb = ism{9b) = 1// (64 in the figure), when the stable and unstable fixed points meet at 
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a single point = bQ^Pb), where the diagonal is tangent to the map (this bifurcation is 
sometimes also referred to as a saddle node bifurcation [2Z])- 

Although in the vicinity of the bifurcation point the behavior is universal, the existence of 
irreversible collisions changes the global dynamics of the map, which in general is not chaotic. 
For example, for angles slightly larger than 6h, the map exhibits the universal intermittent 
behavior associated with tangent bifurcations, with long-lived intervals of quasi-periodic 
motion of period g = 1, as the impact parameter goes across the 60 ~ bl region [27j (see 
ref. [2Hj for a discussion of the corresponding renormalization-group approach in a generic 
tangent bifurcation, and ref. [21] for a more recent perspective on tangent bifurcations in 
the context of Tsallis statistics). 

However, the behavior outside this quasi-periodic region is also periodic due to repeated 
irreversible collisions. A typical trajectory with such intermittent behavior is shown in fig. 
[Sjo. The trajectory is quasi-periodic with g = 1 for 69 ~ K- other hand, outside 

the near critical region, the trajectory is also periodic after two irreversible collisions, with 
period g = 9 as shown in figures |3]d and |4j 

In general, there are two types of irreversible collisions, which we may call positive and 
negative ones, with +1 and —1 as the outgoing impact parameter, respectively, independent 
of the value of the incoming parameter. Therefore, a trajectory becomes periodic whenever a 
second positive or negative collision occurs, which again highlights the commensurate nature 
of the periodic orbits. 

We can still use the Poincare map to determine the phase locking behavior for all forc- 
ing angles. In the q = 9 periodic trajectory shown in fig. |3)d, for example, the particle 
goes through the bottom-top periodic boundary condition only once, which means that its 
asymptotic migration angle will be arctan(l/8), as shown in figure [s] Figure [s] also shows 
the behavior near the critical angle 6c and the entire 6 vs. a curve, which exhibits the typical 
Devil's staircase structure [30]. We also show some of the observed locking angles and their 
ordered structure, in which they form a Farey tree [31] . 

C. Limiting Behavior 

In figure [6] we present the impact parameter map for different magnitudes of the external 
force. In the limit of large external forces the repulsive centers have negligible influence on 
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FIG. 5: Migration angle as a function of the forcing angle (same system as in figs. |3]and|4|. The 
critical forcing angle corresponding to a tangent bifmxation is shown with the dashed line. The 
inset shows the bifurcation region in more detail. The open circle corresponds to the forcing angle 
9 = 12.71°, (3 = 0.55 > g = 9 and directional locking into the [1,8] direction, also discussed in 
figs. |3]and|4] 

the trajectories, and the map tends to an indentity relation, with both 6c —^^ 1 and hj ^ Bq. 
On the other hand, for small forces, that is for forces / ~ 1, all collisions become irreversible, 
in that 6c ^ for / ^ 1. Figure [6] shows that the map of the impact parameter tends to a 
piecewise-continuous map with only two regions: the region of unperturbed trajectories for 
|6o| > 1, and the region of irreversible collisions for |6o| < 1- This is a particularly interesting 
limit, in that it corresponds to impenetrable repulsive cores which could model, for example, 
the transport of a tracer particle through an array of impermeable obstacles in the limit of 
high Peclet numbers, as we investigate in more detail in the next section. 
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FIG. 6: Poincare map of the incoming impact parameter into itself for various relative magnitudes 
of the external force /. The dimensionless lattice parameter is £ = 2.5. 

IV. STOCHASTIC TRANSPORT: HIGH PECLET NUMBER BEHAVIOR 

In this section we investigate the stochastic transport of colloidal particles in the periodic 



landscape discussed in section IIB i. e., we consider diffusive transport in addition to 
the purely convective motion considered in the deterministic case. We will show that, 
at relatively high Peclet numbers, the average motion of the particles exhibits directional 
locking equivalent to that observed in the deterministic case. 

In the presence of diffusive transport, the effective migration angle is given by the angle 
between the average velocity of the particles and the periodic lattice [32] • We first solve the 
Fokker-Planck equation for the probability density associated with the stochastic motion 
of the colloidal particles given by Eq. [l| The Fokker-Planck equation in non-dimensional 
variables reduces to: 

|p(x, t) + f^P{^, t) - 1^ V^P(x, t) = 0, r > 1, (11) 
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FIG. 7: Average migration angle calculated from Eq. 12 for different forcing angles. The periodic 
landscape is described in section II B| The dimensionless parameters are: I = 4.0 and / = 1.0. For 
all cases we computed the steady state probability distribution using a finite element method with 
over 3 x 10^ de grees of freedom and an element size smaller than 10 ^ x R close to the obstacle 
boundary. 

where the Peclet number is given by Pe = FmaxR/D'y, with D the diffusivity of the col- 
loidal particle. The asymptotic distribution of colloidal particles in the unit cell, Poo(x), 



corresponds to the steady state solution of equation [TT] above inside the unit cell and us- 
ing periodic boundary conditions We can then obtain the components of the average 
migration velocity, (Ui), as well as the migration angle a from their ratio [33j, 

(Uy) 



tanfo;) 



dy (-DVPoo(x)) 



dx (/Poo(x)-DVPoo(x)) 



L^O 



(12) 



We obtained the stationary solution for the probability distribution using standard nu- 
merical methods. In figure [7] we present the results for the average migration angle as a 
function of the forcing angle for different Peclet numbers. It is clear that, as the Peclet 
number increases and the convective transport becomes dominant, the relation between the 
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FIG. 8: Asymptotic distribution of particles in the unit cell for different Peclet numbers. The 
dimensionless parameters are: i = 4.0 and / = 1.0. The forcing angle is 9 = arctan0.15 ~ 8.53°. 
The plot at the center shows the concentration profile on a line along the y axis (perpendicular 
to the external force). The dashed line in the plot marks the intersection with the tangent at the 
bottom of the circle. Both lines are shown in the probability plot corresponding to Pe=500. 
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migration and the forcing angle tends to a structure similar to those observed in the deter- 
ministic case, consisting of plateaus and steps. In fact, we show that the average effective 
angle exhibits, in the limit of high Peclet numbers, the same directional locking as that 
predicted in the deterministic case (see figure [?]). 

In figure [8] we present the asymptotic probability distribution in the unit cell for different 
Peclet numbers. We can see that, as the Peclet number increases and convective transport 
becomes dominant, the probability flux is dominated by a stream that leaves the obstacles 
from the bottom of the circle and is parallel to the force {x axis). This is consistent with our 
interpretation in the deterministic case. Note that the solutions presented here correspond 
to the limiting case / = 1.0 and, therefore, there is no penetration into the obstacles in the 
analogous deterministic discussed in section [III C[ 

The plot in figure [8] shows the concentration profile on a line that is perpendicular to 
the forcing direction. It is clear that, as the Peclet number increases, the probability flux 
focuses on a narrow region that streams from the bottom point of the obstacle. In fact, 
the probability maximum in the cross-section plot in figure [8] tends to the point at which 
the tangent to the circle parallel to the force intersects the cross-section line (indicated by 
a dashed line in the plot). 

In the plots corresponding to Pe=100 and 500 we can see the re-entrance of the probability 
stream on the left boundary of the unit cell, as well as its collision with the obstacle. We 
can estimate the width of the probability peak to validate that it corresponds to diffusive 
spreading in the direction perpendicular to the convective motion. For this we follow the 
approximations done in the case of geometric ratchets[3l], in which diffusion along the 
direction of the external force is neglected. In this case, 

(Ay^) = 2DAt~2(^)D=^, (13) 

where £/2 is the distance between the streaming point and the cross-section measurement. 
This estimate agrees well with the numerical results. For example, for Pe=500 the previous 
Eq. predicts Ay ~ 0.09 and the width of the peak is (Xy ~ 0.06. 

Finally, in figure [9] we compare the results of our model to the case of impermeable solid 
obstacles. We discussed in the previous section that, in the limit of / ~ 1 all collisions are 
irreversible {be = for / = 1) and do not penetrate the obstacles. In this case we expect 
the dynamics to be similar to that in the case of solid obstacles. This is only approximate, 
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FIG. 9: Migration angle as a function of the forcing angle. Open symbols correspond to the 
parabolic repulsive centers with non-dimensional lattice spacing I = 4.0, / = 1, and Pe=1000. 
Solid symbols corresponds to numerical simulation in the case of solid, non-permeable obstacles 
for the same lattice spacing and Peclet number. 

given that in the presence of diffusion, there is a non-zero probabiUty of finding the particles 
inside the repulsive regions. Figure |9] shows that, in fact, the migration angles obtained in 
our model are very similar to the migration angles obtained in the case of solid obstacles for 
the same, relatively large, Peclet number. 

V. RELEVANCE TO MICROFLUIDIC DEVICES 

In experimental work, it has been shown the possibility to separate particles using differ- 
ences in the value of the bifurcation angle. In dimensional variables, the bifurcation angle 
is given by Lsin(6'f,) = R/f, which shows its dependence on the properties of the particles. 
Specifically, the properties of the particles could come into play through the force ratio / or 
through the effective size of the repulsive centers R. 
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In the case of optical lattices, both the magnitude of the repulsive force as well as the 
characteristic size of the repulsive centers depend on the particle size [TT]. In fact, our results 
for the bifurcation angle are analogous to the calculation by Pelton et al. [H] for the critical 
angle to escape a single barrier. Our results can then be used to investigate the behavior at 
angles above the bifurcation. 

Finally, we can investigate the transport of suspended particles in a pattern of solid 



obstacles using the limiting case / ~ 1 discussed in III C 



In this case, sin(^fe) = R/L, which depends on the size of the particles through the 
effective size of the obstacles, R = Rq + a, where Rq is the size of the obstacles and a is 
the radius of the particles. In this trivial approximation the bifurcation angle can also be 
obtained by straightforward geometrical considerations. The size dependency implies that 
larger particles will become unlocked from the a = direction at larger angles of the driving 
force, as observed in experiments. 

Consider, for example, the separation of particles of two different sizes, ai = 4/im and 
02 = 6/im, in a lattice of (solid) cylindrical obstacles of radius Rq = 5/im and with a 
lattice constant L = 25^m. Let us assume that the particles are driven at an average 
velocity of f/ ~ lOfim/s, which corresponds to a large value of the particles Peclet number, 
Pe~ 10^. Then, considering only hard-sphere interactions (excluded volume effects) between 
the spheres and the obstacles the corresponding effective size of the obstacles is Ri = 9/xm 
and i?2 = for particles of radius ai and a2, respectively. Similarly, the dimensionless 

length depends on the size of the particles. Specifically, ii = 25/9 and £2 = 25/11, which 
results in different bifurcation angles for the two type of particles. Obi = 21.1° and 6b2 = 26.1° 
for sizes ai and 02, respectively. Therefore, for a driving force oriented at any angle 6 relative 
to the obstacle lattice such that Obi < 9 < 9b2, the particles will separate. In particular, 
the small particles will move in the lattice direction [1,2] with ai = 26.56° and the large 
particles will be locked at 0:2 = 0°. 



VI. SUMMARY 



We have shown that the transport of particles in a periodic lattice of repulsive centers 
exhibits analogous behavior to that observed in microfluidic systems. The simplicity of the 
parabolic repulsive potentials allowed us to integrate the trajectories explicitly and showed 
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the existence of periodic trajectories that are commensurate with the obstacle lattice. We 
also showed that the motion can be determined by means of a Poincare map of the incoming 
impact parameter into itself, which shows that there is a tangent bifurcation at the critical 
forcing angle for which the locking becomes different from a = 0. The entire migration-angle 
vs. forcing-angle curve shows the characteristic Devil's staircase type of structure common 
to phase-locking systems. Finally, we showed that the limiting behavior for impenetrable 
obstacles is equivalent to the high Peclet number limit in the case of transport of particles 
in a periodic pattern of solid obstacles. Therefore, our previous results provide insight into 
the separation problem in the case of periodic potential landscapes as well as in the case of 
periodic patterns of solid obstacles. In fact, we discuss a straightforward application of our 
results to calculating the bifurcation angle in both solid and optical lattices. 

VII. ACKNOWLEDGEMENTS 

We thank R. Hansen for pointing to us the Farey-tree structure underlying the possible 
locking angles. This material is partially based upon work supported by the National Science 
Foundation under Grant No. CBET-0731032. 



[1] K. D. Dorfman and H. Brenner, J. Colloid Interface Sci. 238, 390 (2001), ISSN 0021-9797. 
[2] K. D. Dorfman and H. Brenner, Phys. Rev. E 65, 052103 (2002), ISSN 1063-651X. 
[3] L. R. Huang, E. C. Cox, R. H. Austin, and J. C. Sturm, Science 304, 987 (2004), URL 
|http : //www. sciencemag. org/cgi/ content /abstract /304/5673/987 



[4] J. P. Beech and J. O. Tegenfeldt, Lab chip 8, 657 (2008). 

[5] K. J. Morton, K. Loutherback, D. W. Inglis, O. K. Tsui, J. C. Sturm, S. Y. Chou, and 



R. H. Austin, Proc. Natl. Acad. Sci. pp. 7434-7438 (2008), URL |http : //www . pnas . org/cgi/ 
|content/abstract/0712398105vl| 
[6] K. J. Morton, K. Loutherback, D. W. Inglis, O. K. Tsui, J. C. Sturm, S. Y. Chou, and R. H. 



Austin, Lab chip 8, 1448 (2008), URL |http : //dx . doi . org/10 . 1039/b805614e! 
[7] P. T. Korda, M. B. Taylor, and D. G. Crier, Phys. Rev. Lett. 89, 128301 (2002), ISSN 0031- 
9007. 



21 



[8] K. Ladavac, K. Kasza, and D. Grier, Phys. Rev. E 70 (2004). 

[9] A. Gopinathan and D. G. Grier, Phys. Rev. Lett. 92, 130602 (2004). 

[10] Y. Roichman, V. Wong, and D. G. Grier, Phys. Rev. E 75 (2007). 

[11] M. Pelton, K. Ladavac, and D. Grier, Phys. Rev. E 70 (2004). 

[12] C. Reichhardt, C. J. Olson, and M. B. Hastings, Phys. Rev. Lett. 89 (2002), ISSN 0031-9007. 

[13] C. Reichhardt and F. Nori, Phys. Rev. Lett. 82, 414 (1999), ISSN 0031-9007. 

[14] V. I. Marconi, S. Candia, P. Balenzuela, H. Pastoriza, D. Dommguez, and P. Martinoli, Phys. 
Rev. B 62, 4096 (2000), URL |http : //link . aps . org/abstract/PRB/v62/p4096l 



[15] C. Reichhardt, C. J. O. Reichhardt, and M. B. Hastings, Phys. Rev. E 69, 056115 (2004), 
ISSN 1063-651X. 

[16] C. Reichhardt and C. J. O. Reichhardt, Phys. Rev. E 69, 041405 (2004), ISSN 1063-651X. 
[17] E. Ott, Chaos in Dynamical Systems (Cambridge University Press, 2002), 2nd ed., ISBN 
9780521010849. 

[18] S. H. Strogatz, Nonlinear dynamics and Chaos: with applications to physics, biology, chem- 
istry, and engineering (Addison- Wesley Pub., Reading, Mass., 1994), ISBN 0201543443, URL 
|http : //www. loc . gov/catdir/enhancements/f y083d/93006166-ci.html 



[19] A. Lacasta, J. Sancho, A. Romero, and K. Lindenberg, Phys. Rev. Lett. 94, 160601 (2005). 
[20] J. M. Sancho, M. Khoury, K. Lindenberg, and A. M. Lacasta, J. Phys. Condens. Matter 17, 

S4151 (2005), ISSN 0953-8984. 
[21] A. M. Lacasta, M. Khoury, J. M. Sancho, and K. Lindenberg, Mod. Phys. Lett. B 20, 1427 

(2006), ISSN 0217-9849. 



[22] R. Guantes and S. Miret-Artes, Phys. Rev. E 67 (2003), URL ,http : //link . aps . org/ 
|abstract/PRE/v67/e046212[ 

[23] J. P. Gleeson, J. M. Sancho, A. M. Lacasta, and K. Lindenberg, Phys. Rev. E 73, 041102 
(2006), ISSN 1539-3755. 

[24] D. A. Mc Quarrie, Statistical Mechanics (University Science Books, 2000), 2nd ed., ISBN 
9781891389153. 

[25] T. F. Hannes Risken, The Fokker-Planck Equation: Methods of Solutions and Applications 

(Springer, 1996), 2nd ed., ISBN 9783540615309. 
[26] V. I. Arnold, Ceometrical Methods in the Theory of Ordinary Differential Equations (Springer, 

1996), 2nd ed., ISBN 9780387966496. 

22 



[27] J. P. Eckmann, Rev. Mod. Phys. 53, 643 (1982). 

[28] B. Hu and J. Rudnick, Phys. Rev. Lett. 48, 1645 (1982). 

[29] F. Baldovin, Phys. A 372, 224 (2006). 

[30] P. Bak, Phys. Today 39, 38 (1986). 

[31] H. G. Schuster, Deterministic Chaos (VCH, New York, 1995). 
[32] Z. Li and G. Drazer, Phys. Rev. Lett. 98, 050602 (2007), ISSN 0031-9007. 
[33] H. B. David Edwards, Macrotransport Processes (Butterworth-Heinemann, 1993), ISBN 
9780750693325. 

[34] C. Keller, F. Marquardt, and C. Bruder, Phys. Rev. E 65, 041927 (2002), ISSN 1063-651X. 



23 



